Global Analysis of Fragmentation Functions for Pions and Kaons 

and Their Uncertainties 
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We present new sets of pion and kaon fragmentation functions obtained in NLO combined analyses 
of single-inclusive hadron production in electron-positron annihilation, proton-proton collisions, and 
deep-inelastic lepton-proton scattering with either pions or kaons identified in the final state. At 
variance with all previous fits, the present analyses take into account data where hadrons of different 
electrical charge are identified, which allow to discriminate quark from anti-quark fragmentation 
functions without the need of non trivial fiavor symmetry assumptions. The resulting sets are in 
good agreement with all data analyzed, which cover a much wider kinematical range than in previous 
fits. An extensive use of the Lagrange multiplier technique is made in order to assess the uncertainties 
in the extraction of the fragmentation functions and the synergy from the complementary data sets 
in our global analysis. 
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I. INTRODUCTION AND MOTIVATION 

The hadronization process turns partons produced in 
hard-scattering reactions into the physical, colorless, non 
perturbative hadronic bound states towards any hard 
interaction ultimately develops. Within the standard 
framevirork of leading-twist coUinear QCD [ij, processes 
with an observed hadron in the final-state can be de- 
scribed in terms of perturbative hard scattering cross 
sections and certain non-perturbative but universal func- 
tions: parton distributions, accounting for the partonic 
structure of the hadrons in the initial state just before 
the interaction, and fragmentation functions, encoding 
the details of the subsequent hadronization process . 

These three ingredients are therefore the pillars of the 
perturbative QCD (pQCD) description of hard interac- 
tions; their precise knowledge has been crucial for its 
success in the past, and is imperative for the ongoing 
and future high energy research programs In the 
last few years, the improvement in each of these key ar- 
eas has been remarkable. Higher order QCD calculations 
have been explored and validated for most processes up 
to next-to-leading order (NLO) accuracy, and are cur- 
rently being extended even beyond that point for some 
observables Q. 

The knowledge on parton distributions has become in- 
creasingly precise as a result of two decades of a wide 
variety of high precision measurements, and strenuous 
efforts to update and enlarge periodically the correspond- 
ing QCD analyses [1, Q. State-of-the-art sets of parton 
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densities agree with each other well within the already 
fairly small estimated uncertainties, and provide a pic- 
ture of the proton structure fully consistent with the 
data. For most observables, the differences arising from 
the use of one or another modern set of parton densities 
are negligibly small compared with the uncertainties in 
their measurement or from unknown higher order correc- 
tions. 

Also fragmentation functions have been rapidly evolv- 
ing, following the path of parton densities, however, with- 
out attaining yet the precision of the latter 0, H, H, [l^ ■ 
Most of the information used to determine these distri- 
butions comes essentially from electron-positron annihi- 
lation into charged hadrons. These data have the im- 
portant advantages of being very precise, thanks to high 
statistics measurements from CERN-LEP and at SLAC 
[O, [H, in, Q El mi and clean, in the sense that the 
corresponding cross sections have no dependence on par- 
ton densities. 

In spite of these advantages, electron-positron anni- 
hilation data suffer from many shortcomings. In the 
first place, the data give per se no information on how 
to disentangle quark from anti-quark fragmentation as 
they always refer to the charge sum for a certain hadron 
species, e.g., 7r+ + tt^. The information on how the indi- 
vidual quark fiavors fragment into hadrons depends cru- 
cially on so-called "tagging" techniques and the under- 
lying assumptions implemented in the Monte-Carlo gen- 
erators employed [l4]. In addition to that, at the mass 
scale of the Z-boson, i.e., for the bulk of the electron- 
positron annihilation data, all electro- weak couplings be- 
come roughly equal and thus only flavor singlet combi- 
nations of fragmentation functions can be determined. 

Also the gluon fragmentation is not exceedingly well 
constrained, since the sub-leading NLO corrections for 
electron-positron annihilation are too weak to determine 
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it; quark-gluon mixing in the scale evolution of fragmen- 
tation functions is not enough of a constraint either, due 
to the lack of precise enough data at energy scales away 
from the Z-resonance. The data also become consid- 
erably less accurate and sparse at large hadron energy 
fractions, leading to large uncertainties in that region. 

Fortunately, in the last few years several one-particle 
inclusive measurements coming from both proton-proton 
collisions and deep-inelastic lepton-nucleon scattering 
have matured enough as to yield complementary infor- 
mation on the fra gme ntation process with competing pre- 
cision l2ll . [2^ . These measurements not only 
probe fragmentation in a complementary energy regime, 
but also weigh quite differently contributions of individ- 
ual parton flavors in the hadronization process. 

Including these data in the extraction of fragmentation 
functions, not only increases statistics but, most impor- 
tantly, yields a much more complete picture of the frag- 
mentation process. For example, the complementary in- 
formation allows to relax and test certain rather stringent 
assumptions on the flavor symmetry of fragmentation 
functions usually made. One can also make an indepen- 
dent check of the flavor tagging techniques implemented 
in electron-positron annihilation, and, ultimately, scruti- 
nize the fundamental ideas of QCD factorization and the 
universality of the fragmentation functions. 

In this paper, we perform for the first time such a more 
comprehensive global QCD analysis to obtain new sets of 
pion and kaon fragmentation functions in agreement with 
the wealth of electron-positron annihilation data and also 
electron-proton and proton-proton observables. The lat- 
ter were either impossible to calculate with previous sets 
or just not reproduced by them. As it should be expected, 
our new sets agree with previous extractions 0, M, B 
in aspects of fragmentation functions which are actually 
determined well by electron-positron data but show sig- 
nificant differences otherwise. 

In order to asses the uncertainties in the resulting frag- 
mentation functions, associated with both the uncertain- 
ties in the data and the theoretical estimate of the ob- 
servables, we have made an extensive use of the Lagrange 
multiplier method [2^. We study in detail the profiles 
of the x^-function that quantifies the agreement between 
the set of fragmentation functions, more specifically their 
moments, and the data. This procedure gives a clear pic- 
ture of the relative uncertainty for each quark flavor and 
the gluon, which are found to be about 3% and 10% for 
"favored" (valence) fragmentation functions in pions and 
kaons, respectively, and of the order of 10% and 20% 
for the respective "unfavored" (sea) fragmentation func- 
tions. The method allows also to assess the role and 
interplay of the various data sets in constraining the dif- 
ferent fragmentation functions, illustrating the synergy 
characteristic of a global analysis. 

In the following Section, we briefly summarize the 
QCD framework for fragmentation functions and the dif- 
ferent single-inclusive hadron production processes in- 
cluded in the global analysis. In Sec. Ill, we outline the 



details of the analysis, discussing our choice for the func- 
tional form used to parametrize the fragmentation func- 
tions at the initial scale of the scale evolution, and the 
data sets included in the fit. We also outline the imple- 
mentation of the Mellin transform method for a fast eval- 
uation of the NLO cross sections, and the Lagrange mul- 
tiplier technique for assessing uncertainties. In Sec. IV 
we discuss in great detail our results for both pion and 
for kaon fragmentation functions and their uncertainties. 
We briefly summarize our results in Sec. V. 



II. QCD FRAMEWORK FOR 
FRAGMENTATION FUNCTIONS AND 
SINGLE-INCLUSIVE HADRON PRODUCTION 

The pQCD framework for single-inclusive hadron pro- 
duction in electron-positron annihilation, lepton-nucleon 
deep-inelastic scattering, and hadron-hadron collisions is 
in place for quite some time now, and calculations at 
NLO accuracy are "state-of-the-art" throughout. In each 
case one exploits the factorization theorem [ij, which 
states that the cross section can be decomposed into ap- 
propriate convolutions of perturbatively calculable par- 
tonic hard scattering cross sections and certain combina- 
tions of non-perturbative parton distribution and frag- 
mentation functions. We can restrict ourselves to a brief 
summary of theoretical framework relevant for our global 
QCD analysis of fragmentation functions to set up our 
notation in which we closely follow Refs. and [2^. 



A. Properties of Fragmentation Functions 

A field theoretical definition of fragmentation functions 
in terms of bi-local operators was given in Ref. 
and reads for quarks (up to kinematical pre-factors) 
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and similarly for gluons. V and V' denote the necessary 
gauge- links to render ([1]) gauge-invariant. The simple 
parton model interpretation of (z) as the probability 
for a parton / to produce a hadron H with fraction z 
of its momentum is recovered in light-cone gauge where 
V = V = 1. Since a specific hadron H with light-cone 
momentum is observed in the final-state, a local op- 
erator product expansion (OPE) does not apply. For 
fully inclusive parton densities, the OPE is the basis for 
first principle computations of some of their integer mo- 
ments within "lattice QCD" . Similar calculations cannot 
be pursued for fragmentation functions. 

The scale dependence of the fragmentation functions 
is calculable in pQCD and governed by renormal- 
ization group equations very similar to those for parton 
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densities. For instance, the singlet evolution equation 
schematically reads 

d 



dlnQ2 
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p(T) ^ 



(T) 



(T) 



(2) 



(3) 



(4) 



gg 



is the matrix of the singlet timelike evolution kernels. 

(T) 

The NLO splitting functions have been computed 

in [i^jll^ or can be related to the corresponding spacelike 
kernels by proper analytic continuation |,2S] . 

The range of applicability for fragmentation functions 
as defined above is severely limited to medium-to-large 
values of z. On the one hand, the timelike evolution ker- 
nels in (j4]) develop a strong singular behavior as 2: ^ 0, 
and, on the other hand, the produced hadrons are con- 
sidered to be massless. More specifically, the splitting 
functions Pgq (z) and Pgg {z) have a dominant, large 
logarithmic piece ~ In z/z in their NLO part, which 
ultimately leads to negative fragmentation functions for 
z ^ 1 in the course of the evolution and, perhaps, to 
unphysical, negative cross sections, even if the evolution 
starts with positive distributions at some scale Qo < Q- 
At small z, also finite mass corrections proportional to 
Mn/isz"^) become more and more important. While 
there are ways to resum the singular small-z behavior 
to all orders in as, there is no systematic or unique way 
to correct for finite hadron masses, for instance by intro- 
ducing some "re-scaled" variable z' in SIA. Inseparably 
entwined with mass effects are other power corrections 
or "dynamical higher twists" . 

Anyway, including small- 2; resummations or mass cor- 
rections in one way or the other in the analysis of hadron 
production rates is not compatible with the factoriza- 
tion theorem and the definition of fragmentation func- 
tions outlined above. "Resummed" or "mass corrected" 
fragmentation functions should not be used with fixed or- 
der expressions for, say, the semi-inclusive deep-inelastic 
production of a hadron, eA*" e'HX, discussed in 
Sec. Ill CI Therefore we limit ourselves in our global 
analysis to kinematical regions where mass corrections 
and the influence of the singular small-z behavior of the 
evolution kernels is negligible. It turns out that a cut 
z > ^min = 0.05 (0.1) is Sufficient for data on pion (kaon) 
production. 

Finally, conservation of the momentum of the frag- 
menting parton / in the hadronization process is sum- 
marized by a sum rule stating that 
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H 



f dzzDf{z,Q^)^l, 
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(5) 



i.e., each parton will fragment with 100% probability into 
some hadron H. Equation ([5]) is compatible with the evo- 
lution kernels in the MS scheme, although not for each 
individual contribution dzzDf {z,Q^). Of course, the 
sum rule ([5|) should be dominated, perhaps almost sat- 
urated, by the fragmentation into the lightest hadrons 
such as pions and kaons. The unstable small- 2: behav- 
ior, however, prevents Eq. ([5]) from being a viable con- 
straint in a global analysis. Only truncated moments 
dzzD^^ {z,Q^) are meaningful. 



B. Single-inclusive e^e Annihilation 

The cross sections for the single-inclusive e+e" anni- 
hilation (SIA) into a specific hadron H, 



e e 



(7, Z) 



(6) 



at a center-of-mass system (c.m.s.) energy ^/s and in- 
tegrated over the production angle can be written as 



H 



1 da 
(Ttot dz 



[2F^{z,Q^)+FS{z,Q^)\. (7) 
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The energy Eh of the observed hadron scaled to the beam 
energy Q/2 — is denoted by z = 2pH ■ q/Q^ = 

2Eh/\/s with Q being the momentum of the intermedi- 
ate 7 or Z boson. 



(8) 



is the total cross section for e+e^ hadrons including 
its NLO 0{as) correction and ctq ~ \.'Kc?{Cf')l s. The 
sums in ([7]) and ([5]) run over the n f active quark flavors q, 
and the Cg are the corresponding appropriate electroweak 
charges (see App. A of Ref. [2J] for details). 

To NLO accuracy, the unpolarized "time-like" struc- 
ture functions F^^ and F^ in ([7]) are given by 

2Ff(z,Q^) = Y.eli[Df{z,Q') + D^{z,Q')] 
g I 




(9) 



g 

+Cg^®Z?f] (z,Q2), (10) 

with (8) denoting a standard convolution. The relevant 
NLO coefficient functions C;l'i" in the MS scheme can be 



found in App. A of Ref. [2j 
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We note that the longitudinal structure function Fl 
in Eq. receives its leading nonzero (finite and scheme 
independent) contribution at 0{as). We treat, however, 
the 0{as) expressions in (fTU)) as sub-leading (=NLO) 
in calculations of the total (longitudinal plus transverse) 
cross section ([T]). For predictions of only the longitu- 
dinal cross section at NLO, the C'(a^) corrections [26] 
should be included. However, such measurements are 
not available for identified pions or kaons considered in 
this analysis. 

C. Semi-Inclusive Deep-Inelastic Scattering 



with the NLO (MS) coefficient functions C/^^ [HIMllO, 
p. 

In our global analysis of fragmentation functions we 
will make use of (preliminary) data for charged pion and 
charged kaon multiplicities taken by the HERMES exper- 
iment dl]. The multiphcities {l/NBis)dN" /dzdQ^ are 
defined as the ratio of the semi-inclusive deep-inelastic 
scattering (SIDIS) cross section pT|) in a certain bin of, 
say, and z, to the totally inclusive DIS rate. The par- 
ticular value of this data in the global analysis emerges 
from the sensitivity to individual quark and anti-quark 
flavors in the fragmentation process which is not accessi- 
ble from e+e^ annihilation. 



The cross section for the semi-inclusive deep-inelastic 
production of a hadron, eiV —>■ e' HX, is proportional to 
certain combinations of both the parton distributions of 
the nucleon N and the fragmentation functions for the 
hadron H. It can be written in factorized form in a way 
very similar to the fully inclusive DIS case [2^.[29l.[30l.[3l|: 



da 



H 



dx dy dzff 



y 



y 



2F,"{x,ZH,Q^) 

(11) 



with X and y denoting the usual DIS scaling variables 
{Q^ = sxy), and where \2§, ^3d\ zh = Ph ■ Pn/pn ■ 1 
with an obvious notation of the four-momenta, and with 
—q^ = . Strictly speaking, Eq. (fTTj) and the vari- 
able zji only apply to hadron production in the current 
fragmentation region. This is usually ensured by a cut 
xp > on the Feynman- variable representing the frac- 
tional longitudinal c.m.s. momentum. If necessary, target 
fragmentation could be accounted for by transforming to 
the variable |3l|,[13| zh ^ z = En^i-x) ' energies Eh, 
E]\j defined in the c.m.s. frame of the nucleon and the vir- 
tual photon, and by introducing the so-called "fracture 
functions" [32]. 

The structure functions F^ and in (fTTj) are given 
at NLO by 



2Fl'ix, ZH, Q") = J24\ Q^)D"izH, Q') 



2^ r 



F^ix,ZH,Q^)^ 




qg 



ix,ZH,Q^), (13) 



D. Hadron-Hadron Collisions 

The single-inclusive production of a hadron H at high 
transverse momentum px in hadron-hadron collisions is 
also amenable to QCD perturbation theory. Up to correc- 
tions suppressed by inverse powers of pT, the differential 
cross section can be written in factorized form as 1251 13311 



fb da, 



ab ' 



D 



H 



(14) 



where the sum is over all contributing partonic channels 
a + b —t c + X, with rfa^^ the associated partonic cross 
section, rfijjj^ can be expanded as a power series in the 
strong coupling as and the 0{a^g) NLO corrections are 
available [2^,|33|. As always, the factorized structure ([14]) 
forces one to introduce into the calculation scales of the 
order of the hard scale in the reaction - but not specified 
further by the theory - that separate the short- and long- 
distance contributions. We have suppressed the explicit 
dependence on these renormalization and factorization 
scales in Eq. (|14|). for details, see, e.g., Ref. ^25*]. 

In studies and quantitative analyzes of hadronic cross 
sections, NLO corrections are of particular importance 
and generally indispensable in order to arrive at a firm 
theoretical prediction for ([H]) . Since NLO corrections 
are known to be significant, LO approximations usually 
significantly undershoot the available data. In addition, 
hadronic reactions suffer from much enhanced theoretical 
uncertainties than the reactions described above due to 
the presence of more non-perturbative, scale dependent 
functions. The dependence on the unphysical factoriza- 
tion and renormalization scales can be only controlled 
and quantified at NLO (or beyond). 

As will be discussed below, the special value of 
hadronic cross sections in a global analysis of fragmenta- 
tion functions is their enhanced sensitivity to the gluon 
fragmentation function through the dominance of gg —> 
gX processes for hadrons produced at low-to-medium 
transverse momenta and their sensitivity to fragmenta- 
tion at very high z. Charge separated data for H = tt 
and provide additional information on the flavor sep- 
aration of the . 



± 



5 



III. OUTLINE OF THE ANALYSIS 

In this Section, we outhne the details of our analysis. 
More specifically, we discuss our choice of parametriza- 
tion, the selection of data sets, treatment of experimental 
normalization uncertainties, and how we determine the 
parameters by means of a global minimization. We 
also briefly sketch how we make use of Mellin moments 
to include exact NLO expressions for the cross sections 
(O , (|lip , and in our analysis and how we assess un- 
certainties in the extraction of fragmentation functions 
with the help of the Lagrange multiplier technique. 



A. Parametrization 

All recent analyses of fragmentation functions are 
based exclusively on SIA data 0, H, S and have 
chosen the most simple functional form A^iz"'(l — z)^^ 
to parametrize the at some initial scale /xq for the 
Q^-evolution ^ . The structure of the SIA cross section 

([7|)- (fTU)) allows to extract only information on D'^^^^ 
from data (similarly for kaons). Without assumptions it 
is impossible to distinguish "favored" or "valence" from 
"unfavored" or "sea" fragmentation, for instance, D^^ 
from ZJg where |7r+) = \ud). This is a serious limita- 
tion of all present analyses [3, [l^ , as the obtained 
fragmentation functions cannot be used to compare to a 
wealth of recent data on the production of charged pions 
and kaons in SIDIS or proton-proton collisions [2]| . 
In Ref. a linear suppression factor / = {1 — z) 
was assumed to break this "deadlock" . This was later 
shown to be in fair agreement with charged pion multi- 
plicities in SIDIS from HERMES [3 within a LO com- 
bined analysis of SIA and SIDIS data [111; see also Fig. [3] 
and discussions below. 

In our global analysis we will determine for the first 
time individual fragmentation functions for quark and 
anti-quarks for all flavors as well as gluons from data. 
To accommodate also the experimental information from 
lepton-nucleon and hadron-hadron scattering data, we 
adopt a somewhat more flexible input distribution than 

in [3,[i,i,[i3 

. ^ ^.z"-(l~z)ft[l-f7.(l-^)^-] 
' ^ ' B[2 + a„ ft -fl] + 7»5[2 + a,, ft + ft + 1] ' 

(15) 

where B[a,b] represents the Euler Beta- function and Ni 
is normalized such to represent the contribution of D^^ to 
the sum rule ([5]). A more restrictive initial parametriza- 
tion with 7i = in Eq. (|15[) would introduce artifi- 
cial correlations between the behavior of fragmentation 
functions in different regions of z obscuring also the as- 
sessment of uncertainties. We find that the extra term 
~ (1 — z)*^* in Eq. (jlSp considerably improves the quality 
of the global fit, closely related to the fact that the anal- 
ysis of fragmentation functions is restricted to medium- 



to-large z. Accordingly, additional power terms in z, em- 
phasizing the small z region, have little or no impact on 
the fit and are not pursued further. The initial scale /iq 
for the Q^-evolution is taken to he fig = 1 GeV in our 
analysis. 

Since the initial fragmentation functions (jlSp at scale 
fiQ should not involve more free parameters than can be 
extracted from data, we have to impose, however, cer- 
tain relations upon the individual fragmentation func- 
tions for pions and kaons. We have checked in each case 
that relaxing these assumptions indeed does not signif- 
icantly improve the of the fit of presently available 
data to warrant any additional parameters. In detail, for 
{u, u, d, d} — > tt"*" we impose isospin symmetry for the 
sea fragmentation functions, i.e., 

Df^Df, (16) 

but we allow for slightly different normalizations in the 
q + q sum: 

Df+i-N^ZU- (17) 
For strange quarks it is assumed that 

Df = Df = N'Dl^ (18) 

with N' independent of z. 

It is worth noticing that assuming N = N' — 1 0, [l3| 
in Eqs. P?|) and p^ . respectively, SIA data alone allow 
to distinguish between favored and unfavored fragmenta- 
tion functions in principle. We shall scrutinize the com- 
patibility of these assumptions with SIDIS and hadronic 
scattering data in Sec. IIVFI At any rate, their impact 
on the assessment of uncertainties of fragmentation func- 
tions is highly non trivial. 

For charged kaons we fit D^^^ and Df_^^ independently 
to account for the phenomenological expectation that the 
formation of secondary ss pairs, which is required to form 
a \K^) = I us) from a u but not from an s quark, should 
be suppressed. Indeed, we find from our fit, see Sec. 
IV below, that -D^g > 0^^+^^ in fine with that expec- 
tation. For the unfavored fragmentation the data are 
unable to discriminate between flavors and, consequently, 
we assume that all distributions have the same functional 
form: 

We adopt the functional form (jlSp also for the fragmen- 
tation of heavy charm and bottom quarks into charged 
pions and kaons but setting 7^ = 0. As in 0, [1, we 
assume that Df = and Df = for H = Tr+ , K+ . 
Heavy flavors are included discontinuously as massless 
partons in the evolution ^ above their MS "thresholds" , 
Q = rricfi, with rric^b denoting the mass of the charm and 
bottom quark, respectively. This treatment of heavy fla- 
vors is very much at variance with heavy quark parton 
densities, where very elaborate schemes have been devel- 
oped to properly include mass effects near threshold and 
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to resum large logarithms ~ Inm^^/Q^ for ^ "^c6- 
Only SIA data at ^/s 3> mc.b are sensitive, however, to 
charm and bottom fragmentation in the analysis. Nei- 
ther the charged pion or kaon multiplicities in SIDIS 
nor hadron production data from proton-proton collisions 
at RHIC receive any noticeable contribution from heavy 
quark fragmentation. Therefore the massless approxima- 
tion outlined above, also adopted in 0, H, d, [l3| , appears 
to be sufficient for the time being. However, we note that 
a dynamical, parameter-free generation of the heavy fla- 
vor component to light meson fragmentation functions 
based on NLO matching conditions has been developed 
recently in 35| . This might prove to be a viable alterna- 
tive to the presently adopted framework in the future. 

Thus in total we have to determine 23 (24) parameters 
in the global analysis describing the hadronization of 
quarks and gluons into positively charged pions (kaons). 
Corresponding fragmentation functions into tt^, are 
obtained as usual by charge conjugation and those for 
neutral pions by assuming D'^ — [D^ + ]/2. 

We numerically solve the renormalization group equa- 
tion in NLO 



din /i^ 



Po 2/ 2\ 



f3i 



;al{^?) (20) 



with /3o = 11 — 2ny/3 and /3i = 102— 38n//3 to determine 
the running of the strong coupling. The number of active 
flavors Uf IS increased upon crossing the heavy flavor 
thresholds at /i^ = m\ and = to^ for which we choose 
[3^ TTT-c = 1.43 GeV and nih — 4.3 GeV. We specify Us in 

l("/=4) _ 



the solution of (P(7| for nj = 4 by adopting [36| Aq^j-, 
334 MeV. For our leading order (LO) analysis we use [3^ 

in Eq. ^ 



A^"/-**^ = 220 MeV and, of course, set (3x 



B. Selection of Data Sets 

The parameters describing the fragmentation functions 
for pions and kaons at scale /xq in Eq. are determined 
by a standard minimization for N data points, where 



N 

2=1 



(r. - E,f 

SE? 



(21) 



Ei is the measured value of a given observable, SEi the 
error associated with this measurement, and Ti is the 
corresponding theoretical estimate for a given set of pa- 
rameters in (|15l) . Since the full error correlation matrices 
are not available for most of the data entering the global 
analysis, we take, as usual 0)[l)S[i3|) the statistical and 
systematical errors in quadrature in 6Ei. 

In (PT|) we use charged pion and kaon production data 
in SIA from TPCjIl at = 29 GeV, SLD 0], ALEPH 
[ni, DELPHI [ifTand OPAL [H, all at = Mz- To 



further constrain fragmentation functions through scale 
evolution, we also use data from TASSO (17i] at interme- 
diate c.m.s. energies of ^/s = 33 and 44 GeV, which suffer, 



however, from rather large experimental uncertainties as 
compared to the other SIA data listed above. Other mea- 
surements of SIA [3^1 have too large experimental uncer- 
tainties and hence are not used in our analysis. Because 
of the conceptual problem with fragmentation functions 
at small z outlined in Sec. Ill Al the cut Zmin = 0.05 (0.1) 
is imposed for all pion (kaon) data sets. 

Besides these fully inclusive measurements also "fla- 
vor tagged" SIA results are available, where the quark 
flavor refers to the primary qq pair created by the in- 
termediate photon or Z-boson. ALEPH DELPHI 
p^ . and TPC [l^ provide tagged results distinguish- 
ing between the sum of light u, d, s quarks, charm, and 
bottom events. This information is of particular value 
for the flavor decomposition, as the fully inclusive, "un- 
tagged" data mainly constrain the flavor singlet combi- 
nation (z) on the Z-resonance due overwhelming sta- 
tistical precision of the LEP and SLD data and the fact 
that ~ at Q = ^z- On the downside, flavor tagged 
results can neither be measured directly nor can they be 
unambiguously interpreted and calculated in pQCD. Fla- 
vor enriched samples are unfolded based on Monte-Carlo 
simulations estimating the flavor composition of the data 
sets. For heavy flavor tagged data, a further complica- 
tion arises due to possible contaminations by weak decay 
channels. In our analysis, we obtain the corresponding 
theoretical results Ti by summing in Eqs. (|7|)- (fT0l) only 
over those flavors which are tagged experimentally. At 
NLO this accounts for gluon radiation as well as the pos- 
sibility that not the original (anti-)quark but the radiated 
gluon produces the observed hadron. 

In addition to the flavor tagged results just discussed, 
OPAL [13] has presented fully flavor separated "data" 
in terms of "probabilities" rj^{xp,s) for a quark flavor 
i — q -\- q to produce a "jet" containing the hadron H 
with a momentum fraction z larger than Xp. Needless to 
say, that these results are even more difficult to interpret 
within pQCD beyond the LO and should not be taken 
to literally. To take this into account we assign an up 
to 10% normalization uncertainty to the OPAL tagging 
probabilities 77/^ |1J] in the fit. Nevertheless some tension 
with other data sets remains, in particular for 77^ and 77^, 
as will be discussed below. In our analysis we interpret 
the OPAL results as 



vfixp,Q = Mz 



dz- 



1 da 



H 



cjtot dz 



(22) 



where the subscript i — q denotes that in Eqs. (|7))- (fTm) 
all sums only include the specific quark flavor i. 

To further constrain the fragmentation of different fla- 
vors, as well as too separate favored (valence) and un- 
favored (sea) fragmentation, we include experimental in- 
formation from SIDIS, see Sec. Ill Cl More specifically, we 
make use of (preliminary) charged pion and kaon multi- 
plicities from the HERMES experiment [I8]. These data 
also provide an important consistency check of the flavor 
separation obtained from flavor tagged SIA "data" , as 
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well as of pQCD scale evolution ^ since they refer to 
much lower scales ~ Q = 1 3 GcV <C Mz ■ 

In the minimization we have to account for the fact 
that the SIDIS data are taken in certain bins of z and 

, whereas the theoretical estimates (fTTj) are computed 
for the center of each bin. We estimate the correspond- 
ing uncertainty as the maximal variation of the cross sec- 
tion within each particular bin, see Figs . Bl and [T51 below . 
and add it in quadrature to SEi in Eq. pT|) . In case of 
charged kaon multiplicities, we allow for an additional 
5% uncertainty to account for a possible inadequacy of 
the massless approximation. 

A wealth of new data on single-inclusive hadron pro- 
duction from RHIC experiments [l^, [13, Hi], [22l| have also 
been included in our global analysis. These encompass 
the pt spectrum of neutral pions at central rapidities 
\v\ < 0.35 by PHENIX ^ and at three different forward 
rapidities (r?) = 3.3, 3.8, and 4.0 by STAR [20]. For the 
latter we exclude the most forward bin, (77) = 4.0, from 
the fit as it has large theoretical uncertainties due the 
small Pt values probed. BRAHMS has very recently pub- 
lished Pt spectra for identified charged pions and kaons 
at two values of (forward) rapidities rj = 2.95 and 3.3 [2l| 
of which we use only the former in the fit for similar rea- 
sons as above for STAR. In addition, there are data on 
K% production at central rapidities jryj < 0.5 from STAR 
[23 |. To accommodate the Kg data in the fit, we assume 
that i^g = {K+ + K-)/2 with u ^ K+ and d ^ K+ 
fragmentation functions interchanged. 

For all hadronic data from RHIC, an additional 5% 
error is assigned in quadrature to 5Ei in evaluations of 
in Eq. (PT|) as a rather conservative estimate of the 
theoretical uncertainties related to the choice of the fac- 
torization and renormalization scales in (jl4p . 

To ease possible tensions between certain data sets, we 
allow the data to "float" within the normalization uncer- 
tainties quoted by each experiment. More precisely, in 
addition to the 0(20) parameters describing the frag- 
mentation functions in (|15p , we also fit a set of relative 
normalization factors for each experiment in the mini- 
mization to determine the optimum fit. We note that the 
possibility of normalization uncertainties has been not 
addressed in all previous analysis of SIA data [7..i^.[9l.[l0j. 



C. Mellin Technique 

The integro-differential evolution equations ^ can be 
straightforwardly solved analytically in Mellin n-moment 
space along the lines described, e.g., in Ref. [3§|. The 
Mellin moments of, for instance, the fragmentation func- 
tions Df{z, Q^), are defined as 

Df(n,Q2)= f dzz^-^D^{z,Q^), (23) 
Jo 

and can be expressed in terms of Euler Beta functions for 
our ansatz pS]) at scale /io. The relevant moments of the 



(T\ 

evolution kernels P^^ (n) are given in |40||. The evolved 
fragmentation functions in z-space are re-obtained by an 
inverse Mellin transform given by 



Df{z,Q' 



1 
27ri 



dnz'"Df{n,Q^), 



(24) 



where C„ denotes an appropriately chosen contour in the 
complex n plane. 

The property that numerically very time-consuming 
convolutions in z-space factorize into simple products un- 
der Mellin moments makes them also an ideal tool to 
compute cross sections. Also "plus distributions" , which 
regularize singularities as 2: — > 1, are much easier to han- 
dle. For the SIA cross section ([7]) - ([TU]) the virtue of 
Mellin moments is immediately obvious since they can 
be taken analytically for the hard scattering coefficient 
functions C^'g and can be found in Ref. [40]. We note 
that also the "tagging probabilities" r]f in Eq. (|^^ as ob- 
tained by OPAL JJi] can be straightforwardly computed 
in Mellin moment space 



dn 



1 — xi} 1 da 



H 



atot dn 



l=q 

(25) 

where \i-a denote the Mellin moments of the SIA 

cross section ([7|) - (fT0| for a single flavor i = q. 

The direct use of the SIDIS cross section (fTT1) -(fT3 |) in 
our global analysis would be rather time-consuming and 
awkward though not impossible since the partonic coeffi- 
cient functions Clj^ are still fairly simple. Again, trans- 
forming Eqs. (fTT |) - (fT3|) to Mellin space is much more ap- 
propriate in extensive numerical analyses. As for SIA, the 
Mellin moments of the coefficient functions can be taken 
completely analytically and can be found in [2^, l4l| . 
Since the Cj-^ depend both on x and z, a double Mellin 
transform is required for SIDIS: 



4^'^^(n,m) 



= / dx: 



dz : 



X, z ) 



(26) 

where the dependence on the factorization and renormal- 
ization scales is suppressed in (|26p . Upon combining the 
C^j'^\n,m) with the appropriate n and m moments of 
the evolved parton densities fi{n, Q^) and fragmentation 
functions Dj^ {m,Q^), respectively, the SIDIS cross sec- 
tion is obtained by a numerical fast double inverse Mellin 
transform as described in [41 ] . 

For the much more complex and lengthy partonic hard 
cross sections entering a calculation of hadronic cross sec- 
tions, Eq. p4|) . at NLO accuracy, Mellin moments can no 
longer be computed analytically. Nevertheless it is in the 
analysis of hadron-hadron collision data where the Mellin 
moment technique exhibits its full potential and useful- 
ness [41I, [4^. The crucial, but simple, "trick" in applying 
Mellin moments to Eq. p4p is to express the fragmen- 
tation functions (z) by their Mellin inverses {n) 
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defined in Eq. ([M]) . One subsequently interchanges inte- 
grations and arrives schematically at 



dPH 27ri ^ Jc^ 

X ^ fa fb da^t, 



(27) 



One can now pre-calculate the quantities da^f^{n) = 
J2abfa-'^fb'^ d<^ab ® z'"- in (p7)l , which do not depend on 
the fragmentation functions (n), prior to the fit for 
each contributing combination of single-inclusive subpro- 
cesses producing a certain parton c and in each experi- 
mental bin. We emphasize that in this way all the tedious 
and time-consuming integrations are already dealt with. 

The inverse Mellin transformation, which finally links 
the moments of the fitted fragmentation functions with 
the pre-calculated da'^^{n), of course still needs to be per- 
formed in each step of the fitting procedure. However, 
the integration over n is extremely fast to perform by 
choosing the values for n on the contour C„ in ((27)l sim- 
ply as the supports for a Gaussian integration. The point 
here is that the integrand in n falls off very rapidly as |n| 
increases along the contour for two reasons: first, each 
fragmentation function is expected to fall off at least as 
a power (1 — z)*^, /3 > 1, at large z, which in moment 
space converts into a fall-off of ~ or higher. Sec- 

ondly, we may choose a contour C„ that is bent by an 
angle with respect to the standard vertical direction such 
that for large |n|, {z)~^ decreases exponentially along the 
contour |4l| . This greatly improves the numerical con- 
vergence of the calculation of the (i(T^j(n) in Eq. ([27|) and 
also gives them a rapid fall-off at large arguments. We 
note that for all practical purposes of our global analysis 
between 64 and 92 n moments of da'^f^{n) are sufficient to 
reproduce the cross section ([HI) to an accuracy of much 
better than 1% for all data points used in the fit. 

The crucial asset of the Mellin method is the speed 
at which one can calculate the full hadronic cross sec- 
tion at NLO without approximations, once the moments 
da'^i^{n) have been pre-calculated. To give an example, a 
full NLO computation of all 78 data points from RHIC 
used in the analysis of the pion fragmentation functions 
takes much less than 0.1 second as compared to several 
minutes using Eq. pi)) directly. Since a few thousand 
evaluations of each data point are required in course of 
the minimization, this clearly shows the value of using 
Mellin moments. 



D. Uncertainties: Lagrange Multiplier Technique 

The most difficult but crucial issue to be addressed 
in a global analysis is the estimate of the uncertainties 
in the extraction of the individual fragmentation func- 
tions D!^ . Without a proper assessment of errors, any 
interpretation of the results of the fit or predictions for 



observables based on the fitted Df are incomplete and 
perhaps even misleading. Uncertainties in global analy- 
ses have been thoroughly studied in the context of un- 
polarized parton distributions (PDFs) [23,, j56., ^] , where 
the number and precision of the data available is much 
more significant. A reliable estimate of the errors arising 
from PDFs in predictions for observables related to, e.g., 
new physics or Higgs boson searches at CERN-LHC is of 
utmost importance. 

The possible sources of uncertainties for parton den- 
sities or fragmentation functions can be classified into 
those associated with experimental errors on the data, 
and those related to theoretical or phenomenological as- 
sumptions in the global fitting procedure. The latter 
include, for example, higher order QCD effects in the an- 
alyzed cross sections and their scale dependence, the par- 
ticular choice of the parametric form of the distributions 
at the initial scale, and other model assumptions such as 
flavor and charge conjugation symmetries. Clearly, while 
the first category is usually under control, the second one 
is particularly difficult to quantify. 

Many strategies have been conceived and explored in 
order to assess the uncertainties of PDFs and their prop- 
agation to observables, specially those associated with 
experimental errors in the data. These include the "Hes- 
sian approach" [23] , which assumes that the deviation in 
for the global fit is quadratic in the parameters speci- 
fying the input PDFs away from their optimum fit values. 
Then one propagates these uncertainties of PDFs linearly 
to observables. Alternatively, the "Lagrange multiplier 
method" [11] probes the uncertainty in any observable 
or quantity of interest much more directly. It relates the 
range of variation of one or more physical observables 
dependent upon PDFs to the variation in the x^ used to 
judge the goodness of the fit to data. Specifically, it can 
be implemented by minimizing the function 



(28) 



with respect to the set of parameters {aj} describing 
the PDFs, for fixed values of the Lagrange multipliers 
Xi. Each one of the parameters Aj is related to an ob- 
servable Oi depending on {aj}. The choice = in 
p8| corresponds to the optimum global fit {aj}, for 
which X^ii^j}) = Xo ^-nd Oi{{aj}) = 0°. Minimizing 
^{Xi, {aj}) for Xi deteriorates the quality of the 
fit to data and other values for the observable Oi{{aj}) 
are found from the set of newly fitted parameters {aj}. 
From a scries of global fits for different values of Xi, the 
X^({%}) profile depending on different values of Oi can 
be mapped out. In other words, this tell us by how much 
the fit to data deteriorates if we force the PDFs to yield a 
prediction for an observable different to the one obtained 
with the best fit 0°. 

The value and practical feasibility of the Lagrange mul- 
tiplier technique has been demonstrated not only for the 
hig hly sophisticated global analyses of unpolarized PDFs 
(23| . but also in case of polarized PDFs 44]. In Sec. lIVFi 
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we will show that the same holds for the analysis of frag- 
mentation fmictions. Here, the limitations due to the 
available data are in some sense similar to those we en- 
counter for polarized PDFs: the bulk of the data (SIA 
and spin-dependent DIS in case of fragmentation func- 
tions and polarized PDFs, respectively) neither deter- 
mine the giuon well nor allow for a reliable flavor separa- 
tion. In both types of analysis, SIDIS and hadronic data 
provide invaluable constraints on the parameter space 
describing the input densities. 

In an ideal situation, where every source of uncertainty 
is well understood and fully accounted for, all profiles, 
including those for the parameters {aj} of the fit, would 
be parabolic, and the 1-ct uncertainty for any observable 
would correspond to a change in by one unit, i.e., 
Ax^ = 1. This is, of course, rarely the case, and in order 
to account for missing correlated experimental errors or 
theoretical uncertainties in global analysis it is customary 
to consider instead of Ax^ = 1 a 2-^5% variation in as 
a more conservative estimate of the range of uncertainty 

In addition to the possibility of assessing the uncer- 
tainties of parameters aj or observables Oi{{aj}), the 
Lagrange multiplier method allows to elucidate the role 
of each subset of data included in the fit in constraining a 
certain quantity. One just needs to determine the shape 
and variation of the partial contribution of a par- 
ticular subset n of data to the total x^ as the observable 
changes depending on the Lagrange multipliers. When a 
given subset of data can by itself determine, say, a given 
observable Oi{{aj}), the profile of Ax^ w.r.t. Oi{{aj}) 
is expected to be roughly parabolic, with a minimum 
close to the "preferred" value determined by the op- 
timum global fit. However, when a given subset of data 
does not fully constrain the observable Oi{{aj}), its pro- 
file w.r.t. Oi{{aj}) is either fiat or increases (decreases) 
monotonically without minimum in the range of variation 
of the observable. In general, constraints on Oi{{aj}) in a 
global analysis result from the subtle interplay of several 
subsets of data. The combination of the different partial 
contributions to x^, even of those that by themselves do 
not show a minimum, define the final x^ profile and the 
best fit value 0°, thereby highlighting the complemen- 
tary nature of a global analysis. For our fragmentation 
functions this will be illustrated in detail in Sec. IIVFI 



IV. RESULTS 

In this Section we discuss in detail the results of our 
global analysis of fragmentation functions for pions and 
kaons. We present the parameters of the optimum fits 
describing the at the input scale /zq, compare 

to the data used in the analysis, and give x^ values for 
each individual set of data used. Detailed comparisons 
are made with the results obtained in the analyses of SIA 
data in Refs. and @, in the following labeled as KRE 
and AKK, respectively. Even though we are mainly inter- 



ested in a precise extraction of fragmentation functions at 
NLO accuracy, we also briefly present corresponding re- 
sults of a global analysis performed at LO approximation. 
The significantly better x^ of the NLO sets highlights the 
importance of the NLO corrections and the limitations of 
a LO analysis. Nevertheless, our LO sets should be used 
in calculations of observables where NLO corrections are 
not available, or in event generators limited to LO accu- 
racy. 



A. NLO analysis of pion fragmentation functions 

From a first glance at Figures[l][51 one immediately no- 
tices the remarkable agreement between our new NLO fit 
and data. Experimental results for inclusive hadron pro- 
duction in SIA and proton-proton collisions span several 
orders of magnitude, and the energy scale of the different 
processes ranges from 1 GeV to the mass of the Z-boson. 
This strongly supports the underlying theoretical frame- 
work outlined in Sec. [TTl in particular the fundamental 
notions of factorization and universality for fragmenta- 
tion functions. 

Existing sets of NLO pion fragmentation functions 
0, 0] also give a nice overall description of the SIA 
data included in these analyses, as indicated in Figs.[T][31 
They fail, however, to satisfactorily reproduce charged 
pion production data obtained in SIDIS and in proton- 
proton collisions. Figs. S] and [HI respectively. In addition, 
estimates for neutral pion production rates in proton- 
proton collisions based on KRE [T^ or AKK 'y\ frag- 
mentation function differ substantially as can be seen in 
Figs. [5] and [71 On the contrary, our new set of NLO 
fragmentation functions gives, for the first time, a nice 
global description of hadron production data in electron- 
positron, lepton-nucleon, and hadron-hadron scattering, 
which constitutes a significant and necessary improve- 
ment. 

The most significant difference between our NLO 
global analysis and previous extractions of Df in 0, d, H, 
[ig | is the fact that we can now determine most aspects of 
the fragmentation functions from data rather than being 
forced to make assumptions due to the insufficient infor- 
mation contained in the SIA data alone. We find that, 
in particular, the extra freedom regarding flavor symme- 
try (or the lack thereof) as introduced in Eqs. P7)) and 
p8)) allows us to reproduce the charged pion data shown 
in Figs, m and [Sj In Table [J we give the set of param- 
eters specifying the optimum flt of pion fragmentation 
functions at NLO accuracy in Eq. (fT5)) at our input scale 
/io = 1 GeV for the light flavors and the gluon, and at 
fiQ ~ nic = 1.43 GeV and /io = rrih — 4.3 GeV, for charm 
and bottom fragmentation, respectively. As can be in- 
ferred from there, the outcome of the global analysis de- 
viates from the symmetry assumptions 0, [l3| N = 1 and 
TV' = 1 in Eqs.dni) and ^ by more than 10% and 20%, 
respectively. 

Another crucial asset of our analysis is the enhanced 
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FIG. 1: L.h.s. comparison of our NLO results for {l/atot)da^ /dz according to Eq. ([7]) with the data sets for inclusive pion 
production in SIA used in the fit, see Tab. HIl R.h.s. "(data-theory) /theory" for our NLO results for each of the data sets. 
Also shown are the results obtained with the KRE [7| and AKK Q parameterizations, dashed and dotted lines, respectively. 




FIG. 2: Same as in Fig. [T] but now for light quark ("uds") tagged cross sections. 
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TABLE I: Parameters describing the NLO fragmentation 

functions for positively cfiarged pions, DJ {z, fio), in Eq. (llSp 
at the input scale no — 1 GeV. Inputs for the charm and bot- 
tom fragmentation functions refer to fio = mc ~ 1.43 GeV 
and ^0 = rni, = 4.3 GeV, respectively. 



flavor i 






ft 


7i 


5^ 


u + u 


0.345 


-0.015 


1.20 


11.06 


4.23 


d + d 


0.380 


-0.015 


1.20 


11.06 


4.23 


u — d 


0.115 


0.520 


3.27 


16.26 


8.46 


s + s 


0.190 


0.520 


3.27 


16.26 


8.46 


c + c 


0.271 


-0.905 


3.23 


0.00 


0.00 


b + b 


0.501 


-1.305 


5.67 


0.00 


0.00 


g 


0.279 


0.899 


1.57 


20.00 


4.91 



flexibility of the initial light quark and gluon fragmen- 
tation functions as a function of z in Eq. (|15p as com- 
pared to the standard three parameter form used so far 
01 Hi S • This is not only indispensable to accommo- 
date SIDIS and hadron-hadron scattering data but even 
somewhat improves the quality of the fit to the SIA data. 
Indeed, upon closer examination of Figs. [T][3l in partic- 
ular the "(data-theory) /theory" insets for each data set 
on the right-hand side (r.h.s.) of each plot, one finds a 
slightly improved overall agreement with data as com- 
pared to the, still excellent, one for KRE and AKK; see, 
for instance, the TPC or ALEPH data in Fig. [TJ 

Noticeable is also that our fit follows the trend of the 
data even below the z values included in the analysis. As 
for KRE [7], all data below Zmin — 0.05 are not taken 
into account in the minimization to ensure that the 
possible impact of small-z resummations or hadron mass 
effects, see Sec. Ill Al is negligible. In contrast, the agree- 
ment with data rapidly deteriorates for AKK imme- 
diately below z ~ 0.1, from which they choose to exclude 
data from the fit. This might be linked with the less 
flexible functional form for the fragmentation functions. 

In Tab. |TT] we list all data sets included in our global 
analysis as discussed in Sec. IIII Bl and give the individual 
values for each set. We note that quoted values are 
based only on fitted data points, i.e., z > 0.05 for SIA, 
and include normalization uncertainties determined for 
each experiment in the fit. Allowing for relative normal- 
izations in a global analysis within the range quoted by 
each experiment is a common tool to ease possible ten- 
sions between certain data sets. Indeed we find that the 
global fit considerably improves after taking into account 
normalization "shifts" . 

In spite of the nice "visual" agreement between the fit 
and data found in Figs.[T][8l the total of 843.7 units in 
Tab. [n] appears to be fairly large in view of the roughly 
400 data points fitted. For the SIA data, the large 
can be pinpointed to only very few data points. Due to 
extremely high precision of the data on the Z-resonance, 
any deviation between data and theory is strongly pe- 
nalized in the evaluation and results in an overall 



TABLE II: Data used in the NLO global analysis of pion 
fragmentation functions, the individual values for each set, 
the fitted normalizations, and the total of the fit. 
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o Lag 
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"c tag" 


1.10 


5 


38.3 




"6 tag" 


1.10 


5 


73.0 


HERMES [18] 


^+ 


1.03 


32 


67.4 






1.03 


32 


120.8 


PHENIX [19] 




1.09 


23 


76.4 


STAR [22] 


{ri) = 3.3 


1.05 


4 


3.4 




tt", {-q) = 3.7 


1.05 


5 


9.8 


BRAHMS [21] 


7r+, (r;) = 2.95 


1.0 


18 


28.2 




TT^, (J7) = 2.95 


1.0 


18 


43.0 


TOTAL: 






392 


843.7 



per degree of freedom which is rather large. This is a 
common "characteristic" of all extractions of fragmen- 
tation function made so far 0, H, H, • We also wish 
to point out, that there is a tension between the behav- 
ior of the DELPHI data at large z and those of all the 
other data sets at Q = Mz which cannot be resolved by 
the fit, see Figs. [T]|3l The in general larger values of 
the heavy flavor, in particular bottom quark, tagged SIA 
cross sections in Fig. [3] might be related to some extent 
to contaminations from weak decays. 

We will further scrutinize the quality of the fit to 
the SIDIS and hadronic data in the following. As illus- 
trated in Fig. m the agreement between the (preliminary) 
charged pion multiplicities in SIDIS from the HERMES 
experiment [isj and the results of our fit is remarkably 
good. The theoretical estimates for the multiplicities are 
computed using PDFs from Ref. [3^ as input in Eq. (|12p . 
although no significant differences are found using other 
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FIG. 3: Same as in Fig. [T] but now for charm and bottom quark tagged cross sections. 




1 Q^[GeV] 10 1 Q^[GeV] lO 



FIG. 4: L.h.s. comparison of our NLO results for charged pion muhiplicities in SIDIS, {l/Nnis)dN^ /dzdQ^ , with prehminary 
HERMES data Also shown are the results obtained with the KRE 7] parametrization. R.h.s. "(data-theory)/theory" 

for our NLO results, open and full circles denote tt"*" and vr" multiplicities, respectively. The shaded bands indicate estimates 
of theoretical uncertainties due to finite bin-size effects (see text). 
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FIG. 5: upper panels comparison of our NLO results for single- inclusive charged pion production pp tt^X at rapidity 
T] = 2.95 (solid lines) with BRAHMS data i2l3 using fif — fir = Pt- Also shown are the results obtained with the KRE [3] 
parametrization (dashed lines). The shaded bands indicate theoretical uncertainties when all scales are varied in the range 
Pt /2 < fif = Hr < '2pT- lower panels "(data-theory) /theory" for our NLO results. 



modern sets, e.g., [y]. 

The fit not only reproduces accurately the normaliza- 
tion in each z-bin for both 7r+ and tt", but also the 
"scaling violations" in Q^, which are rather large at the 
low scales involved in the experiment. Using the KRE 
fragmentation functions instead, also nicely reproduces 
the data in the z bins up to 0.6 for the tt"*" multiplic- 
ities but considerably overshoots all 7r~ data, indicat- 
ing that the assumed favored/unfavored separation [7] 
ZJg /DJJ = (1 — z) is not accurate enough. We note 
that AKK refrains from imposing any assumptions in 
their fit, which are beyond what can be obtained from 
SIA data. Hence, the AKK (and also KKP [§]) sets can- 
not be used in theoretical calculations whenever the ex- 
periment discriminates between positively or negatively 
charged pions (or kaons). This, of course, seriously re- 
stricts the potential applications and testability of these 
sets. 

Again, as can be inferred from Tab. [TTl the individual 
from SIDIS is fairly large, in particular for the 7r~ mul- 
tiplicities. As in SIA, these data are based on a sample 
with high statistics and hence the accuracy of the data 
is very good. As before, any deviations between theory 
and data are severely punished by a large contribution 
to x^- It is mainly the bin at the largest z which makes 
all the trouble and contributes most to the total ■ The 
optimum fit already "negotiates" the best compromise in 



describing these data. As mentioned in Sec. IIII Bl we in- 
clude an estimate of the theoretical uncertainties due to 
finite bin-size effects in the minimization. These are 
indicated by the shaded bands on the r.h.s. of Fig.[4l We 
also note, that we do not include the first data point in 
each z bin in the fit as the scale Q almost coincides with 
the already low input scale of our fit. In general, one may 
wonder about possible contaminations from, say, higher 
twists at the low values accessible by HERMES. How- 
ever, we do not find any indications that the SIDIS data 
are incompatible with the other data sets in our global 
analysis. 

In Figs. [5]-[7]we compare the results of our fit to recent 
data from proton-proton collisions at ^/s = 200 GeV at 
RHIC. A general characteristics of all hadronic data is a 
large theoretical uncertainty associated with the choice 
of the arbitrary factorization and renormalization scales 
/i/ and fir, respectively, in Eq. p4p . Although largely 
reduced when going from the LO to the NLO approxi- 
mation for the pp — > ttX cross section, as demonstrated 
in (25l . [3^ , theoretical errors remain much more sizable 
than experimental errors. For our analysis we choose the 
transverse momentum of the produced pion as the hard 
scale, i.e., fij — fir = Pt, which yields a very good overall 
description at both central and forward rapidities. The 
shaded bands in Figs. [5] - [7] indicate theoretical uncer- 
tainties when all scales are varied simultaneously in the 
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FIG. 6: Same as in Fig. [5] but now for single-inclusive neu- 
tral pion production pp -k"X at mid-rapidities < 0.35 
measured by PHENIX Hi]. 



10 



10 

-2 

10 

1 

0.5 


-0.5 
-1 



3 71 

[^ib/GeV^] 




■ THIS FIT 

■ KRE 
AKK 
scale uncertainty 

-I 1 1 1 1— 



(data - theory )/theory 







Pt [GeV] 



FIG. 7: Same as in Fig. |6] but now for STAR data W\ at 
forward rapidities {rj) = 3.3 and (77) = 3.8, solid and open 
squares, respectively. 
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FIG. 8: Comparison of the OPAL "tagging probabilities" [IJ] 
for charged pions, rj^ , as a function of the minimum Xp, see 
Eq. (|22[) . with our NLO results (solid lines). Also shown are 
the results obtained with the KRE and AKK 9,] 

parame- 

terizations, dashed and dotted lines, respectively. 

RHIC data in the global analysis is manifold: data at 
central rapidities 77 ~ and not too large px are strongly 
dominated by gg —^ gX scattering and hence constrain 
the gluon fragmentation Dg. At forward rapidities 77 ~ 3 
the mixture between quark and gluon fragmentation is 
roughly equal. It turns out that for both central and for- 
ward rapidities, the fragmentation occurs at fairly large 
average (z) > 0.5, see, e.g.. Fig. 4 in Ref. where 
information from SIA and SIDIS is sparse. As in case of 
SIDIS, the relevant hard scale of the process, Q — 0{pt), 
is much smaller than in SIA, thereby allowing to exploit 
evolution effects to further constrain the fragmentation 
functions. 

The charged separated pion data obtained by 
BRAHMS very recently and shown in Fig. [5l nicely 
back up the separation of favored and unfavored fragmen- 
tation functions obtained from the SIDIS data discussed 
above. Another important feature of the BRAHMS and 
the RHIC data in general, is the failure of the KRE set 
to reproduce them. As can be seen in Figs. [5] -[71 using 
the KRE fragmentation one considerably undershoots all 
RHIC data. Only by pushing the factorization scales to 
the extreme this could be remedied to some extent. This 
observation has usually been taken as an indication of an 
inadequately small gluon fragmentation function in the 
KRE set at intermediate-to-large z and scales of a few 
GeV. The fact that the agreement with the PHENIX data 
in Fig. [6] is much better at large pt when quark fragmen- 
tation becomes more important, supports this picture. 

The recent AKK set (as well as the preceding KKP 
analysis Q) are characterized by a much larger gluon 
fragmentation function than in KRE and, consequently, 
leads to a good description of the PHENIX data and, to 
a lesser extent, also of the STAR data. The latter may 
suggest the need for an even larger gluon fragmentation 
function. As for SIDIS, the KKP or AKK sets cannot 
be used to compare to the charge separated BRAHMS 
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FIG. 9: Transverse momentum distributions for neutral pions 
obtained by HI [4^ in deep-inelastic e^p collisions compared 
to a NLO prediction using our new set of pion fragmentation 
functions. Also shown are the results obtained with the KRE 
[3 and KKP ^ parameterizations. 



forward and reliable interpretation in pQCD. The fairly 
large values we obtain for yyj, and 77 J are therefore 
not alarming and perhaps even expected. Nevertheless, 
it is encouraging that the general trend of the data is 
roughly reproduced by our combined global fit. We have 
also checked, that the outcome of our analysis does not 
change, if we excluded the OPAL results from out fit. 

Our newly obtained NLO pion fragmentation functions 
are best tested in predictions for cross sections not in- 
cluded in the global analysis. As a final cross check we 
therefore compare in Fig. [5] measurements by the HI col- 
laboration [46| of forward neutral pions in deep-inelastic 
positron-proton collisions at -/s ~ 300 GeV with NLO 
predictions based on our new set of fragmentation func- 
tions. The pions are required to be produced within a 
small angle 6^^ E [5°, 25°] from the direction of the pro- 
ton beam in the laboratory frame, with an energy frac- 
tion = Et^/Ep > 0.01 and transverse momenta in 
the range 2.5 < pr < 15 GeV. The observable has been 
shown to be crucially dependent on NLO contributions 
associated to the gluon fragmentation function [47| . The 
nice agreement between the data and the NLO results 
based on our new set of pion fragmentation functions is 
reassuring. 



data. Since the KRE, KKP, and AKK sets are based 
on the analysis of roughly the same SIA data, the huge 
difference between the obtained gluon distributions only 
demonstrates again that a global analysis is imperative 
in obtaining reliable fragmentation functions. We will 
discuss the individual fragmentation functions and their 
uncertainties in some detail in Sees. IIVEI and lIVFi re- 
spectively. We already note here, that, surprisingly, the 
gluon fragmentation function obtained in our combined 
fit turns out to be smaller than in AKK at intermediate z, 
0-2 ^ z < 0.5, and only becomes larger at z > 0.6. This 
is indicative of the complex interplay of the information 
provided by the different data sets in a global analysis. 

The remaining data set used to constrain the pion frag- 
mentation functions are the OPAL "tagging probabili- 
ties" 14j as defined in Eq. (P^ . As can be inferred from 
Fig. [51 all sets reproduce nicely the "data" for up and 
down flavors, which indicates that they are already essen- 
tially fixed by SIA data. For strangeness fragmentation, 
however, the sets differ considerably. In the KRE anal- 
ysis 0, -Dg is completely fixed by assuming iV' = 1 in 
Eq. (flSl) , which overshoots significantly the OPAL result 
for 77 J. For heavy fiavors, in particular for charm, the 
agreement is less favorable for all sets. The results for 
77^' from OPAL are clearly at odds with the tagged SIA 
data from TPC and SLD shown in Fig. [3] However, as 
explained in Sec. IIII Bi the interpretation of the OPAL 
results beyond the LO should be taken with a grain of 
salt. So is not surprising to find some discrepancy in 
our fit where, contrary to AKK the flavor separation 
comes not only from the OPAL results but also from 
the interplay of other data sets, which have a straight- 



B. LO analysis of pion fragmentation functions 

For completeness, we have also performed a global 
analysis of the same set of data given in Tab. mi where aU 
observables, as, and the scale evolution of the fragmen- 
tation functions are computed at LO accuracy. We use 
the same parametrization (|15p and fltting procedure as 
in the NLO case and outlined in Sec. lIIII The parameters 
of the optimum LO fit are given in Tab. IIIII 

As it can be immediately seen from Tab. llVi the qual- 
ity of the LO fit is significantly worse than the NLO one, 
resulting in a 25% increase in the total x^- Although all 
individual observables show an increase in x^, it is some- 
what more noticeable for proton-proton collision data. 
This is a sensible and expected result as the NLO correc- 
tions are known to be fairly large and important [2^, . 
To make up for the smaller LO partonic scattering cross 



TABLE III: As in Tab.Ubut now describing the LO fragmen- 
tation functions for positively charged pions, {z,jj,o). 



flavor i 






ft 


7i 


5, 


u + u 


0.367 


-0.228 


1.20 


5.29 


4.51 


d + d 


0.404 


-0.228 


1.20 


5.29 


4.51 


u — d 


0.117 


0.123 


2.19 


7.80 


6.80 


s + s 


0.197 


0.123 


2.19 


7.80 


6.80 


c+c 


0.256 


-0.310 


4.89 


0.00 


0.00 


b + b 


0.469 


-1.108 


6.45 


0.00 


0.00 


9 


0.493 


1.179 


2.83 


-1.00 


6.76 
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TABLE IV: Same as in Tab. |TT]but now at LO accuracy. 



experiment 


data 


rel. norm, data points 




type 


in fit fitted 



22.7 
1.9 
5.6 
7.3 
48.1 
21.5 
20.9 
21.3 
9.3 
34.5 
64.4 
45.9 
30.6 
51.9 
20.7 
9.3 
7.5 
66.7 
36.9 
88.8 



HERMES fl8] 


7V+ 




1.03 


32 


53.6 




7r~ 




1.03 


32 


153.9 


PHENIX [19] 


^« 




1.09 


23 


82.2 


STAR [22| tt" 


{v) = 


3.3 


0.95 


4 


15.5 




{v) = 


3.7 


0.95 


5 


11.7 


BRAHMS [21J 7r+, 


{v) = 


2.95 


1.0 


18 


46.3 




iv)^ 


2.95 


1.0 


18 


77.7 


TOTAL: 








392 


1056.8 



sections relevant for RHIC data, the most striking differ- 
ence between the LO and NLO fragmentation functions 
is found for gluons, while the moments for the quark fla- 
vors remain rather stable; cf. Tabs. |T] and IIIII In spite 
of the larger x^, for consistency, our LO sets should be 
used for rough estimates of observables where NLO cor- 
rections are not yet available, or in event generators based 
on matrix elements at LO accuracy. Because of the lim- 
ited usefulness of the LO set, we refrain from going into 
any further details here. 



C. NLO analysis of kaon fragmentation functions 

Our NLO fits to single inclusive kaon production data 
as compared to those for pions, reflect the sensible dif- 
ference in quality of both data sets. Even the most pre- 
cise kaon production data in SIA exhibit experimental 
uncertainties typically twice as large as those found for 
pions. The potentially problematic low-z region is ex- 



TABLE V: Data used in the NLO global analysis of kaon 
fragmentation functions, the individual values for each set, 
the fitted normalizations, and the total of the fit. 



6xp 6riiii6iit 


delta- 


rel. norm. 


data points 


X 




type 


in fi + 

m nt 


nttea 




IFO llol 


incl. 




1 o 


y.o 




incl. 


U.yoo 


lo 


^ A A 

14.4 




uas tag 


U.yoo 


iU 


^ A A 

14.4 




c tag 


n OQQ 

U.yoo 




17 9 

1 / .z 




tag 


n ACQ 

U.yoo 


1 A 


lO.Z 


A T TP'PTT ri 1 1 


incl. 


U.y / 


lo 


1 9 Q 

iz.o 




incl. 


1 n 


1 9 


1 n 

l.U 




'"''uds tag" 


l.U 


1 9 


9 Q 
Z.O 




U LicLg 


1 

l.U 


1 9 




OPAT, [14-1 


u tag 


1 in 


o 


Kj.O 




"d tag" 


1.10 


5 


9.9 




"s tag" 


1.10 


5 


36.8 




"c tag" 


1.10 


5 


44.9 




"6 tag" 


1.10 


5 


18.6 


HERMES [18] 


K+ 


1.03 


24 


15.0 




K- 


1.03 


24 


79.3 


STAR [20] 


Kl 


0.95 


14 


40.0 


BRAHMS [21] K 


^ = 2.95 


1.0 


18 


28.8 


K 


-, = 2.95 


1.0 


18 


21.5 


TOTAL: 






232 


394.1 



TABLE VI: Parameters describing the NLO fragmentation 

functions for positively charged kaons, (z,fio), at the in- 
put scale fiQ = 1 GeV. Inputs for the charm and bottom 
fragmentation functions refer to fio = "^c = 1.43 GeV and 
/io = mi, — 4.3 GeV, respectively. 



flavor i 


Ni 




ft 


7i 


5^ 


u + u 


0.058 


0.705 


1.20 


15.00 


6.02 


s + s 


0.343 


-0.065 


1.20 


4.36 


3.73 


d + d 


0.016 


1.108 


10.00 


10.00 


3.28 


u — s 


0.008 


1.108 


10.00 


10.00 


3.28 


c+c 


0.196 


0.102 


4.56 


0.00 


0.00 


b + b 


0.139 


-0.584 


7.42 


0.00 


0.00 


9 


0.017 


5.055 


1.20 


0.00 


0.00 



pected to set in earlier due to the larger kaon mass. To 
be on the safe side, we raise the cut from z,„i„ ~ 0.05 
to Zmin = 0-1. Clearly, one must expect much less well 
constrained fragmentation functions for kaons. Other- 
wise the minimization proceeds along the same lines 
as for the analysis of the pion data, and the results are 
summarized in Tab.lVl As demonstrated in Figs. [TUl- [TBI 
the overall agreement of our NLO fit and data is again 
remarkably good. The set of parameters specifying the 
obtained kaon fragmentation functions at NLO accuracy 



TPC [15] 


incl. 


0.94 


17 




"lids tag" 


0.94 


9 




"c tag" 


0.94 


9 




"& tag" 


0.94 


9 


TASSO [17] 


incl. (34 GeV) 


0.94 


11 




incl. (44 GeV) 


0.94 


7 


SLD [16] 


incl. 


1.007 


28 




"uds tag" 


1.007 


17 




"c tag" 


1.007 


17 




"b tag" 


1.007 


17 


ALEPH [11] 


incl. 


0.97 


22 


DELPHI [12] 


incl. 


1.0 


17 




"uds tag" 


1.0 


17 




"& tag" 


1.0 


17 


OPAL [13, U} 


incl. 


1.0 


21 




"u tag" 


1.10 


5 




"d tag" 


1.10 


5 




"s tag" 


1.10 


5 
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1.10 


5 




"b tag" 


1.10 


5 
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FIG. 10: As in Fig. [T] but now for inclusive kaon production in SIA. 
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FIG. 11: Same as in Fig. [10] but now for light quark ("uds") 
tagged cross sections. The upper panel shows the comparison 
of the "uds" tagged cross sections with data and the lower 
two panels show "(data-theory) /theory" . 

can be found in Tab. IVII 



Starting with the SIA data in Figs. \W\ - [H 

our fit 

shows again, thanks to the more flexible z-dependence 
of the initial distributions pS]) and the less stringent fla- 
vor symmetry assumptions, a slightly better agreement 
with SIA data than previous analyses 0,0- Also, below 
Zmin = 0.1 the fit still follows the trend of the data, which 
ensures that finite mass effects or the unstable small-z 
behavior of the scale evolution are still negligible. In 
terms of x^, none of the SIA data sets, even the ones 
for heavy flavor tagged cross sections in Fig. [T^ poses 
any problems. The reduced statistical accuracy of the 
kaon data compared to the pion data is immediately ob- 
vious and also reflects itself in larger fluctuations in the 
"(data-theory) /theory" comparisons in Figs. [TUl- [T^ 

The biggest asset of our global analysis is again that 
not only SIA data but, for the flrst time, also SIDIS 
and RHIC data are nicely reproduced as is demonstrated 
in Figs. [T2] - [m It is worth noticing that, at variance 
with what happens for pions, the HERMES SIDIS data, 
Fig. ll3i rule out completely the flavor separation assumed 
in Ref. [7,] for kaon fragmentation functions. The predic- 
tion based on the KRE set overshoots the data by a factor 
of two, whereas our global flt shows much better agree- 
ment. This is in particular true for the multiplici- 
ties, which in SIDIS receive the dominant contribution 
related to the (large) up quark parton density in the pro- 
ton. The production of K~ is linked predominantly to 
strange quarks and anti-up quarks in the proton, both 
of which are much less abundant than up quarks for the 
relatively large momentum fractions x relevant for HER- 
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FIG. 12: Same as in Fig. [10] but now for charm and bottom quark tagged cross sections. 



MES. Here our fit reproduces the magnitude of the 
multiphcities in each z bin, but not the slope. In fact, 
the related is the by far largest contribution to the to- 
tal of the fit, cf. Tab. [V] The theoretical uncertainty 
introduced by the size of the bins in the computation of 
the cross section is also significantly larger than in the 
case of pions. For negatively charged kaons it is further 
amplified by the rapid growth of sea quark densities in 
the proton with Q^. Taking also into account that the 
strange and, to some extent, also the anti-up quark par- 
ton densities are not too well constraint, the agreement 
of the new fit is rather encouraging. 

Our new set of kaon fragmentation functions also yields 
the best description of the BRAHMS [2ll| and STAR 
data t22i]. Figs. [HI and [T5l respectively, which both probe 
the large z regime, z > 0.5, poorly mapped out by SIA 
data. As for pions, the KRE set undershoots all RHIC 
kaon data, whereas AKK overshoots the STAR A'g data. 
As before, the AKK (and KKP) set does not allow to 
compute the charge separated and K~ yields by 
BRAHMS. Again, the discrepancy with the RHIC data 
can be traced back to the behavior of the gluon fragmen- 
tation function at large z and will be discussed further in 
Sec. IIVEI Notice again the very large theoretical scale 
uncertainties for the RHIC data. 

The last data set entering the global analysis of kaon 
fragmentation functions are the OPAL tagging probabil- 
ities ryf [l3| • In view of the already discussed conceptual 
problems with these data, the overall agreement with the 
OPAL results in Fig. [12] is reasonable. The lack of agree- 



ment, most noticeable for charm, again suggests some 
degree of inconsistency with other data sets included in 
the fit. 



D. LO analysis of kaon fragmentation functions 

The LO global analysis of kaon production data yields 
significantly larger values of ■ The parameters describ- 
ing the optimum LO kaon fragmentation functions are 
given in Tab. IVIl] while contributions are in Tab. lVllT] 

The overall increase with respect to the NLO fit is 
about 30%, even larger than the one found in the case 
of pions, and with the most noticeable differences in the 
partial contributions stemming from RHIC and SIDIS 
data. 



TABLE VII: As in Tab.[Vl]but now for the LO fragmentation 
functions for positively charged kaons. 



flavor i 


N, 


ai 


ft 


7i 


5^ 


u-\-u 


0.054 


1.018 


1.20 


15.00 


6.04 


s + s 


0.361 


0.733 


1.20 


20.00 


5.28 


u — s 


0.005 


1.322 


10.00 


10.00 


3.67 


d-f d 


0.010 


1.322 


10.00 


10.00 


3.67 


c + c 


0.214 


0.239 


4.27 


0.00 


0.00 


h + b 


0.147 


-0.464 


7.37 


0.00 


0.00 


9 


0.036 


5.282 


1.20 


0.00 


0.00 
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FIG. 13: Same as in Fig. |4] but now for charged kaons. 




FIG. 14: Same as in Fig. [S] but now for charged kaons. 
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FIG. 15: Same as in Fig. H but now for STAR data '.22] for 
K% production at central rapidities < 0.5. 
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FIG. 16: Same as in Fig. [8] but now for charged kaons. 
E. Fragmentation Functions 



In this Section we shall present an overall description 
of the different fragmentation functions obtained 



in the global fit [48l]and perform a comparison with the 
KRE [7] and AKK % NLO sets based only on SIA data. 

The upper panels of Fig. [T7| show the fragmenta- 
tion functions for positively charged pions at the scale 
q2 = 10 GeV^. As expected, the sum u + u dominates 
over the unfavored distributions u and s. At large values 
of z, there is an important contribution from the valence 
u fragmentation, while at small z the sea distribution 
dominates and u + u ~ 2u. In the same limit, it can 
be observed that the s fragmentation function turns out 



TABLE VIII: 


Same as in Tab. FVlbut now at LO accuracy. 






data 


rel. norm. 


data points 


2 

X 






type 


in fit 


fitted 








incl. 


0.94 


12 


12.7 


TA'^^D [171 


incl. (34 GeV) 


0.94 


4 


2.6 






incl. 


0.983 


18 


18.1 






"uds tag" 


0.983 


10 


21.9 






"c tag" 


0.983 


10 


18.4 






"6 tag" 


0.983 


10 


15.0 


AT.F.PH [1 11 




incl. 


0.97 


13 


14.0 


DELPHI fl2l 




incl. 


1.0 


12 


1.2 






"uds tag" 


1.0 


12 


2.6 






"6 tag" 


1 A 
i.U 


1 


A 1 
4.1 
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"u tag" 
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i.iU 
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1 n 
iU.D 






"s tag" 


1 in 


p: 



oz.y 






"c tag" 


1 1 n 


c; 


00. ± 






"6 tag" 


1 1 n 
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to be smaller than the corresponding u sea distribution, 
as anticipated in Sec. IIV Al when discussing the value of 
the relevant parameter N' . As can also be noticed from 
Fig. I17[ charm and gluon fragmentation are quite sizable 
and comparable to the one of the light quarks at small z. 
This is actually a general feature of heavy quark fragmen- 
tation, opposite to the behavior of the usually less rele- 
vant heavy quark parton distributions. At this scale, the 
bottom channel has not opened yet, but the correspond- 
ing distribution can be observed in Fig. 1181 where the 
same functions are plotted at a higher scale = Mf . 
As expected, heavy quark and gluon fragmentation are 
rather suppressed at larger values of z. 

In the middle and lower panels of Figs. [TTlandfTSl we 
compare our set of fragmentation functions to those from 
KRE and AKK, respectively. The largest differences ap- 
pear for the unfavored quark and gluon distributions and, 
usually, both at large z and near Zmin below which frag- 
mentation functions cannot be used. Notice that for 
AKK, is not available for comparison and that their 
analysis is limited to z > 0.1 rather than z > Zmin — 0.05. 
Since AKK tends to overestimate the SIA cross-section 
outside the fitted region, i.e., below z = 0.1, any dis- 
agreement there is not surprising. 

While the discrepancy with KRE for the light quark 
distributions are reasonably moderate, there happens to 
be a rather large difference at the level of the strange 
fragmentation with AKK. The origin of this disagreement 
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FIG. 17: upper panels individual fragmentation functions for positively charged pions zDf {z,Q^) at = lOGeV'^ for 
i — u + u, u, s + s, g, and c + c. middle panels ratios of our fragmentation functions to the ones of KRE [3]. The dotted 
line indicates the ratio for singlet combination of fragmentation functions zD^ . lower panels ratios of our fragmentation 
functions to the ones of AKK % note that D'^ is not available in the AKK analysis. 




FIG. 18: Same as in Fig. [T7]but now for Q = Mz- Also shown is the 6 + 6 fragmentation function. 
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FIG. 19: Same as in Fig. [17] but now for positively charged kaons K'^ . 
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can be easily understood: SIDIS data, not included in the 
AKK fit, require a smaller fragmentation for u (as can 
be seen in the same plot) and d quarks. Since the singlet 
combination, D-^, is already well constraint by SIA data, 
as can also be seen in Figs. [T71 and [T51 that automatically 
requires an increase in the D'^ fragmentation function. 

The most pronounced differences appear at very large 
values of z for both heavy flavor and gluon fragmenta- 
tion functions. In the heavy quark case, both SIDIS and 
RHIC data are insensitive and SIA data suffer from very 
large experimental uncertainties for z > 0.4 and even 
conflicting results. Therefore, in that kinematical regime 
those densities cannot be determined well, and large dis- 
crepancies between different analyses are expected. For 
the gluon fragmentation, the main source of information 
at z > 0.5 stems from PHENIX data at mid and STAR 
and BRAHMS data at forward rapidity. These are in- 
cluded only in the analysis presented here. BRAHMS 
and STAR data access the largest z values of all data 
sets in our global analysis. As was observed in Figs. \E\ 
- [3 the KRE sets underestimates all experimental data 
(for II f ~ fir = Pt), while the AKK set agrees with 
the PHENIX data but tends to underestimate STAR. 
This explains the differences in Figs. [T71 and \TE\ between 
our gluon fragmentation function and those of KRE and 
AKK. RHIC data tend to favor a significantly larger ^ 
at large z. 

Finally, the comparison between the different distribu- 
tions and ratios in Figs. [T7] and [TH] also highlights the 
importance of the evolution. Even though the dis- 
tributions show a much stronger rise at small z when 
the scale is increased [i^, one also finds a much better 
agreement between the different sets at = M|, where 
most of the very precise SIA data are obtained. In other 
words, the evolution downwards from Mz, were the dis- 
tributions agree best, to scales relevant for RHIC and 
SIDIS data, exacerbate the differences between them. 

Figures [12] and [201 provide the same information and 
comparisons as in the previous ones but for positively 
charged kaons. As expected, the dominant fragmenta- 
tion in the light quark sector corresponds to the strange 
distribution. While heavy quark densities are as large as 
those for pions, the gluon fragmentation turns out to be 
much less sizable, even though it is still larger than those 
from KRE and AKK at large z in order to fit the proton- 
proton data from STAR and, in particular, at forward 
rapidities from BRAHMS. The comparison in the light 
quark sector shows many similarities with the pion case, 
but here the discrepancies are more noticeable with the 
KRE set instead. The global fit requires a smaller contri- 
bution from u quark fragmentation, mostly from SIDIS 
data (distributions like KRE overestimate SIDIS data, 
see Fig. I13p resulting in an increase in the strange sector 
as the singlet is again constrained by SIA data. 



F. Uncertainties 

In order to give a clear and comprehensive picture of 
the typical uncertainties characteristic of the fragmenta- 
tion functions obtained in the global fits, in the present 
Section we apply the Lagrange multiplier technique in- 
troduced in Sec. IIIIDI 

Rather than focusing on the uncertainties of the pa- 
rameters in Eq. (jlSp determining the fragmentation func- 
tions at the initial scale and choosing a particular incre- 
ment Ax^ to judge the quality of the fit, we find it much 
more enlightening to analyze the range of variation of 
other relevant features of the fragmentation functions, 
with a more apparent physical meaning, and take these as 
the characteristic uncertainties of the fit. Notice that the 
range of variation of the fitted parameters are strongly 
correlated; the impact of any of them on the behavior of 
the distributions, or on a given observable, is determined 
also by the values taken by the whole set of parameters 
through the evolution equations. Of course, in order to 
get a precise estimate of the uncertainty in a given ob- 
servable computed with the set, the range of variation 
of that particular observable as a function of Ax^ has 
to be evaluated. As explained in Sec. IIII Dl the result 
takes into account the complex correlations between the 
parameters, implies no assumptions on the profile of x^, 
and allows to consider different Ax^- 

In Figures [m and [2H we show, as an example, the range 
of variation of the truncated second moments of the frag- 
mentation functions 

r,f{xp,Q^)= f zDf{z,Q^)dz, (29) 

for Xp — 0.2 and (3 = 5 GeV, around the values obtained 
for them in the best fit to data, 77/^, against the corre- 
sponding increase in x^- In the lowest order, the second 
moments represent the energy fraction of the parent par- 
ton of fiavor i taken by the hadron H . The truncated 
moments discard the low-z contributions, which are not 
constrained by data and therefore only determined by 
I0W-2; extrapolation. 

As it can be seen in the upper left panel of Fig. [2T1 
the truncated moment niZ.+in associated with -DJJ^^^, is 
the constrained best, with a range of variation of less 
than 3% around the value computed with the best fit, 
assuming a very conservative increase in x^ by 15 units, 
i.e., Ax^ = 15. This comparatively stringent restriction 
comes from the fact that all the observables accounted for 
in the fit have a strong dependence on the correspond- 
ing D^^j^— fragmentation function. Moving to the next 

panel, we find that for the unfavored fragmenta- 
tion function, the truncated moment is less, but still well 
constrained within a 5% range for a similar Ax^- While 
this distribution cannot be determined by SIA data with- 
out fiavor symmetry assumptions, it is constrained by 
SIDIS and proton-proton collisions involving low x con- 
tributions, where the u{x) PDF is large, or at high x 
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FIG. 21: Profiles of for the NLO pion fragmentation fit as a function of the truncated second moments rfl {xp = 0.2, ( 
5 GeV) for different flavors. The moments are normalized to the value r;7o they take in the best fit to data. 



through its aher ego . 

The fragmentation functions for strange quarks have 
much larger uncertainties, and their moment can vary 
by more that 10%, as shown in the upper right panel 
of Fig. EH In our fit, not only D^^ ^ D^'^ but their 
respective uncertainties are found to be quite different 
as well. Notice that these uncertainties in the unfavored 
distributions cannot justify the differences with KRE and 
AKK sets found in the previous Section. 

The moment for the gluon fragmentation function, in 
the lower left panel of Fig. [211 is restricted to vary by less 
than 10%, with the constraint mainly stemming from the 
evolution from the initial scale to the scales relevant for 
each measurement, rather than from a direct contribu- 
tion to a particular cross section. Such kind of contri- 
butions are certainly present in RHIC data, but only in 
narrow intervals of z, so they cannot fix the truncated 
moment for Xp = 0.2 by themselves. For heavy quarks, 
flavor tagged data dominate the fit, and the more pre- 
cise b-tagged data lead to better constrained moments 
for Df^ than for D^l^. 

b+b 

The uncertainties in kaon fragmentation functions. 
Fig. [551 are typically twice as large as those for pions. 



with 10% variations for the total fragmentation func- 
tions containing a favored fragmentation function and 
rather poorly constrained unfavored fragmentation func- 
tions — ~ . The profiles show much 
more significant deviations from the parabolic behavior 
than in the case of pions, and the upper bound for the 
moment of the gluon distribution is much less defined. 

Next we further illuminate the role of the different data 
sets in setting the constraints on the fragmentation func- 
tions. In Figs. [23l and [24l we show the partial contribu- 
tions Axf of a data set i to the increase of from its 
minimum value defined by the best fit against the varia- 
tion of some selected truncated moments. 

In the upper left panel of Fig. [23l it can be noticed 
that both SIA (thick solid line) and SIDIS (thick dashed 
dotted line) data define their own minima for the u -\-u 
moment, slightly to the right and to the left of the best fit 
value, respectively, but within the above mentioned con- 
servative uncertainty range of Ay^ = 15. The other data 
sets fail to develop well defined minima (at least within 
the shown range), and the final result is a compromise 
between all of them. 

Notice that the failure to define a minimum by a given 
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FIG. 22: The same as in Fig. 1211 but now for kaon fragmentation at NLO. 



subset of data may not necessarily imply a weak depen- 
dence of the data on the particular flavor that fragments. 
It may happen because the data cover a limited range in 
z, while the observable we chose picks up contributions 
over a much wider range in z. The seemingly contra- 
dictory "preferences" of two data sets can follow from 
sensitivity to complementary regions in z. 

The truncated moment for the gluon-to-pion fragmen- 
tation function, shown in the upper right panel of Fig. [231 
happens to be well constrained again by SIA and SIDIS 
data, with their respective minima very close to the best 
fit result. For rjl^-g^ in the lower right panel, neither set 
shows a minimum and the final result is a compromise. 
Notice that, here, the most significant contributions are 
those coming from SIA and the OPAL tagging probabil- 
ities. 

Partial contributions Ax| to the truncated moments 
for kaon fragmentation functions. Fig. [24l have less def- 
inite preferences, leading to much larger uncertainties. 
Starting with the moment for s -I- s in the upper left 
panel, we find that SIDIS data now fail to define a min- 
imum. This is mainly due to the limited strangeness 
content in the proton. However, SIDIS contributes with 
OPAL tagged data to balance the preference coming from 



SIA data, which is slightly larger than the best fit value. 
Gluon-to-kaon fragmentation, shown in the upper right 
panel, is mostly constrained by STAR and SIDIS data. 
The u^u moment, receiving contributions from the sup- 
pressed D^' and the doubly suppressed , is nei- 
ther well constrained by SIA nor by SIDIS data, however, 
OPAL and STAR data help to improve the situation. 

Upon the completion of our analysis, a new determi- 
nation of fragmentation functions, again based only on 
SIA data but with a careful assessment of uncertainties, 
was published [l0|. There, it is shown that the large dif- 
ferences found between the most widely used sets (KRE, 
KKP, and AKK) are related to the large z behavior of 
the fragmentation functions, where the uncertainties are 
found to be most significant. Indeed, we can see that 
this is the case for sets based only in SIA data, but the 
situation is considerably improved in a global fit. In the 
left panel of Fig. [21] we show, as an example, the par- 
tial contributions Axf from different sets of data to the 
truncated moment for -D^iji^j, but now taking only large- 
z contributions, Xp = 0.5 in Eq. ([29)1 . Here it can be 
seen that SIA data, as expected, fail to define a mini- 
mum, but the complementary information coming from 
the other sets, dominated by PHENIX and SIDIS, de- 
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to the NLO global analysis of pion fragmentation functions 
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FIG. 24: The same as in Fig. [Ml but for NLO kaon fits. 

fine a clear minimum, witii a well a constrained range of 
variation, again showing the power and importance of a 
global analysis. 

Another interesting difference with the analysis of 
Ref. is that while they accomplish a good fit to SIA 
data under the assumption iV' = 1 in Eq. (|18p . our global 
fit prefers N' ~ 0.83. In order to understand this differ- 
ence, the right panel of Fig.l^shows the partial contribu- 
tion to for different values of N' normalized to A^' for 
our best fit. Although the fit clearly prefers A^' ~ 0.83, 
the uncertainty is large, and for a very conservative Ax^ 
may be taken as marginally consistent with A^' = 1. The 
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FIG. 25: The same as in Fig. [23l but for <+^(a;p = 0.5, Q = 
5GeV) (left panel) and N' (right panel). 

preference for A^' values smaller than one in our global fit 
is driven by the OPAL tagging probabilities, not included 
in [13. 

Finally, we note that projected measurements of 
hadron production in SIA by the BELLE and/or the 
BaBar experiments at low c.m.s. energies would open up 
the possibility for studies of scaling violations with un- 
precedented precision. Such data should help to further 
constrain the fragmentation functions and significantly 
reduce the present uncertainties. 



V. CONCLUSIONS 

We have demonstrated the feasibility of a NLO com- 
bined QCD analysis of single inclusive hadron produc- 
tion data for pions and kaons, coming from electron- 
positron annihilation, deep-inelastic lepton-nucleon scat- 
tering, and proton-proton collisions, collected over a wide 
kinematic range. 

At variance with previous fits based only on electron- 
positron annihilation data, the present analysis includes 
complementary information from other experiments that 
reduce significantly the uncertainties of the resulting 
fragmentation functions. 

In the case of pion fragmentation functions, we find 
that the new SIDIS data provided by the HERMES ex- 
periment, effectively constrain the separation between fa- 
vored and unfavored distributions, a separation that was 
either not implemented in previous sets or it was based 
on certain assumptions. The most recent RHIC results 
provide stringent constraints on the gluon fragmenta- 
tion function and, in general, on the large z behavior 
of the other distributions. For kaons, the new data mod- 
ify significantly the up-to-now standard picture provided 
by previous analyses. Specifically, SIDIS data rule out 
the flavor separation scheme hitherto implemented, while 
RHIC data lead to a new gluon fragmentation function, 
and thus scale dependences, significantly different over 
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the whole range of z. 

The implementation of the minimization in our 
global analysis is numerically fast and efficient and can be 
straightforwardly expanded to any future set of hadron 
production data. With the help of the Mellin moment 
technique, the entire analysis was consistently performed 
at NLO accuracy without resorting to often used approx- 
imations for NLO hard scattering cross sections. 

The success of the global analysis performed here, for 
the first time including observables other than single in- 
clusive annihilation, stands for an explicit check of factor- 
ization, universality, and the perturbative QCD frame- 
work for the description of the corresponding processes, 
providing at the same time much more precise constraints 



on the fragmentation functions. Proton-proton collision 
data and that coming from SIDIS offer a crucial piece 
of information that cannot be disregarded and will be 
increasingly accessible in the near future. 
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